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Abstract 

• It has been established that the entangled polymer dynamics can be reasonably described 

^ ■ by single chain models such as tube and slip-link models. Although the entanglement effect is 

a result of the hard-core interaction between chains, linkage between the single chain models 

■ and the real multi-chain system has not been established yet. In this study, we propose a 
I multi-chain slip-spring model where bead-spring chains are dispersed in space and connected 

I . by slip-springs inspired by the single chain slip-spring model [A. E. Likhtman, Macromolecules 

■ 38, 6128 (2005)]. In this model the entanglement effect is replaced by the slip-springs, not 
^ , by the hard-core interaction between beads so that this model is located in the niche between 
O ■ conventional multi-chain simulations and single chain models. The set of state variables are 
^ i ' the position of beads and the connectivity (indices) of the slip-springs between beads. The 

dynamics of the system is described by the time evolution equation and stochastic transition 
dynamics for these variables. We propose a simple model which is based on the well-defined 
^ I total free-energy and the detailed balance condition. The free energy in our model contains 

00 ■ a repulsive interaction between beads, which compensate the attractive interaction artificially 

generated by the slip-springs. The explicit expression of linear relaxation modulus is also 
derived by the linear response theory. We also propose a possible numerical scheme to perform 
simulations. Simulations reproduced expected bead number dependence in transitional regime 
' between Rouse and entangled dynamics for the chain structure, the central bead difi'usion, and 

the linear relaxation modulus. 
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■ 1 Introduction 



It has been rather established that the entangled polymer dynamics can be reasonably described 
by single chain models where the effect of entanglement is replaced by dynamical constraints such 
as tubes or slip-links [IHl]. For instance, several single chain models [SHTT] have been proposed 
to reproduce viscoelasticity, diffusion, dielectric relaxation, etc, by taking account the relevant 
relaxation mechanisms such as the reptation [TJ[2] , contour length fluctuation [T21[T3] and thermal 
and convective constraint releases PTIITS] . 
Despite the successes of these single chain models, the linkage between the single chain models 
and the real situation with multi-chains has not been clarified yet. There have been lots of attempts 
to extract the parameters for single chain model from multi-chain molecular simulations where the 
entanglement effect is naturally taken into account by the hard-core (excluded volume) interaction 
[T1H22]- Specifically, the primitive path analysis [25H?7] is a realization of the original idea of 
entanglement, and thus, it is promising not only to extract the parameters, also to clarify the 
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microscopic picture (or definition) of the entanglement. Even if some parameters in single chain 
models are extracted by the primitive path analysis, there are still some fundamental difficulties. 
In most cases, the mean-field type single chain description is assumed rather a priori, and the 
assumptions employed in the single chain models have not been fully justified yet. For instance, 
mapping of the cross-correlation between different chains p5H5^ ) to the single chain models has not 
been clarified yet. Thus a model which needs fewer assumptions and based on realistic molecular 
picture is demanding. 

A possible approach in the niche between multi-chain and single chain pictures is the multi-chain 
model where the entanglement is not described by the hard-core interaction but by a coarse-grained 
manner similar to the single-chain models . Actually we have developed such a model called 

the primitive chain network (PCN) model (also referred as NAPLES code) where phantom chains 
are bundled by slip- links to form a network in 3D space [34]. The PCN simulations have been 
performed for several systems such as linear polymers |34 1l37fl39] . symmetric and asymmetric 
stars [38l|40], comb branch polymer [41], polymer blends [38], and copolymers [42], and it has 
been showed that the PCN simulations can reproduce various rhcological properties reasonably. 
Moreover, even under large and fast deformations [151 - H^ the PCN simulations show reasonable 
consistency with experiments. To locate the PCN model between multi-chain and single chain 
models, we have reported a comparison with single chain models on bidisperse blends |50| and 
a comparison with molecular simulations on the network statistics |51j . However, the total free 
energy of the system in the PCN model is not well-defined [52], because the equations describing 
dynamics are not based on the free energy nor the detail balance. (This is because the PCN 
dynamics was modeled rather empirically.) As a result, comparison of some static properties of 
the PCN model with the other models is essentially difficult. 

In this study, we newly propose another multi-chain model based on the entanglement picture, 
which is the multi-chain slip-spring model inspired by the single chain slip-spring model proposed by 
Likhtman several years ago [TU]. Differently from the PCN model, we define the total free energy 
for the new model, and we employ dynamics model (a time evolution equation and stochastic 
processes) which satisfies the detailed balance condition. Thus our model reproduces the thermal 
equilibrium which is characterized by the free energy. In the present paper, we report all the 
equations to construct the model, and also a numerical scheme to implement a simulation code. 
Then some preliminary results for such as the chain dimension, the linear viscoelasticity and the 
center-of-mass diffusion, obtained by the simulations, are reported. 

2 Model 
2.1 Overview 

Figure [1] shows schematic view of the model employed in this study. We consider a network 
composed of bead-spring chains connected by slip-springs in a volume V. We describe the number 
of chains in the system as M, and the number of beads in a chain as N. The beads are connected 
along the chain backbone by linear entropic springs characterized by the average bond size 6, and 
chains essentially behave as ideal, Rouse type chains. Apart from the chain connectivity, some of 
the beads are connected to the other beads by the slip-springs to mimic the entanglement between 
chains. Following the previous works jl01lll| . we treat the degrees of freedom of slip-spring as 
state variables, which obey the Maxwell-Boltzmann type statistics in equilibrium. That is, we 
assume that the equilibrium probability distribution of the state variables is described by the 
Boltzmann weight with the (effective) free energy. The ends of slip-springs move (slip) along the 
chain to reproduce the chain slippage via entanglements. Instead of the continuous sliding dynamics 
proposed in the original single-chain model jlOj . we introduce the discrete hopping dynamics for 
the end of slip-spring between neighboring beads jllj . The slip-springs are stochastically destroyed 
when one of the ends is at the chain ends, and they are also stochastically constructed at the chain 
ends with a certain probability. For the spring force for slip-springs, we employ the linear entropic 
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spring force (The strength of this entropic spring is characterized by the slip-spring parameter 
A''s.) It has been ahcady pointed that shp-hnks or shp-springs effectively give attractive interaction 
between polymers and thus the statistical properties of polymer chains are affected by slip-springs 
[53] . We introduce a repulsive interaction between beads to compensate this attractive interaction 
induced by the slip-springs and to recover the ideal chain statistics. The number density of the 
slip-spring, 0, is another parameter to control the strength of the entanglement effect, which is 
related to the entanglement density. 

The state variables of the system are the position of beads, {Ri^k} (where i and k are the 
indices for chain and bead position on the chain, respectively) , the total number of slip-springs in 
the system Z, and the connection matrix of slip-springs {Sa} (where a is the index for slip-spring 
and the definition of Sa is given later). 



2.2 Equilibrium probability distribution 

First we consider the equilibrium statistical properties. We assume that the chains in our multi- 
chain slip-spring system obey the ideal chain statistics. (There is no interaction between beads 
except the entropic linear springs.) Namely, the probability distribution for the polymer chain 
conformations {Ri^k} is given by 
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Here the subscript "eq" represents the equilibrium quantity. V is the volume of the system, b is 
the bead size, fc^ is the Boltzmann constant, and T is the temperature. 

We consider to put slip-springs in the system, preserving the ideal chain statistics mentioned 
above. (As we mentioned, we assume that the system state, including the slip-springs, is characterized 
by the free energy.) For simplicity, we do not introduce restriction for the slip-spring configurations 
(connectivity). For example, we allow multiple slip-springs to share the same bead, or allow slip- 
springs to connect two beads on the same chain. Because both ends of a slip-spring are attached 
to beads, we need four indices to specify the state of a slip-spring. Hereafter, we describe the fc-th 
bead on the i-th chain as (i, fc), and the state of the a-th slip-spring as Sa = {Sa,i, 5*0,2, Sa,3, Sa,4:) 
where the ends of the slip-spring are located at the bead {Sa,i, Sa.2) a-nd the bead {Sa,3, Sa,4:)- 

We assume that there is no specific interaction between slip-springs. (This assumption is the 
same as one employed in single chain models, where slip-springs behave as one dimensional ideal 
gas pUllllj .) The number of slip-springs is not constant and it is controlled by the effective chemical 
potential for slip-springs [54]. Because slip-springs are statistically independent each other, for a 
given polymer conformation, the probability distribution of the slip-spring state is given by 



P,cii{Sa},Z\{R,,k}) = - 



1 



1 



mR^,k}) z\ 



exp 



E 



3(i2, 



■S„,l,S'a,2 



R 



vZ 



(2) 



where P(X|y) represents the conditional probability of X under a given Z is the total number 
of slip-springs in the system, and N g is the parameter related to the spring constant of slip-spring, 
and V is the effective chemical potential for slip-springs. S({i?i _fc}) is the grand partition function 
of the slip-spring defined as 
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Here } taken for all possible slip-spring indices. Eq ((3)) can be calculated as follows. 
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Substituting eq (|4]) into eq ([2|), we obtain 
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Eqs ([T]) and ((S]) give the equilibrium distribution function of the fuU set of state variables as 
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The effective free energy corresponds to eq ^ is given by 
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The free energy ([7]) consists of several contributions. The first and second terms in the right hand 
side of eq ([7]) are the elastic energies of polymer chains and slip-springs, respectively. The third 
term in the right hand side of eq ([7]) represents the repulsive interaction which compensates the 
attractive interaction caused by slip-springs |53| . The repulsive potential is a soft-core Gaussian 
type and similar to the Flory-Krigbaum potential [55]. However, it should be emphasized that this 
repulsive interaction is not introduced to reproduce the excluded volume effect and the overlapping 
among the chains, but to compensate the artificial attraction by slip-springs. Indeed, the Gaussian 
form comes from the harmonic potential of slip-springs. This repulsive interaction acts not only 
for the beads connected to slip-springs but also for the free beads. Using the effective free energy 
given by ([7]), eq ([6]) can be rewritten as 
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can be defined as 



The equilibrium statistical average, which we describe as (. . . 
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In our model, the number of slip-springs is controlled by the effective chemical potential u. 
However, in practice, the effective chemical potential is not convenient nor intuitive. Thus, instead 
of V, we utilize the average number density of slip-springs, (j) ^ input parameter. (In the 
thermodynamic limit where M, V ^ oo with fixed po, the average slip-spring number {Z)cq — 4>V 
becomes an extensive variable while 4> is an intensive variable.) The equilibrium average number 
density of slip-springs can be calculated from the equilibrium probability distribution given by ([6]) . 
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To evaluate 0, the pair-correlation function of beads is required. This pair correlation function 
can be calculated straightforwardly, since there is no correlation between different chains. Thus 
we have 

M{M ~ 1)7V2 
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Substituting eq pT|) into eq ((TU]), we have 
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Here, po = MN/V is the average number density of beads. Eq (|12p indicates that the average slip- 
spring number depends not only on the effective chemical potential but also on the bead density 
Po and the slip-spring intensity Ns . This is different from the single chain slip-spring model where 
the slip-spring density is depends only on the chemical potential v [TT] . From eq (|12p , the effective 
chemical potential j/ corresponding to a given (f> is obtained as 
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To calculate rheological properties, we need expression of the stress tensor of the system. In 
this work, we employ the following definition for the stress tensor according to the stress-optical 
law. 

1 



V 
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Here 1 is unit tensor. In equilibrium, eq (|14p reduces to the following form. 



)oq - - — ksTl 



(14) 



(15) 



5 



This is the stress tensor of ideal gas, of which number density is M/V . The definition of the stress 
tensor in a shp-spring type model is not trivial. It is possible to include the contributions from 
the slip-springs (contributions from the elastic force and the repulsive force). Fortunately, most 
of rheological properties seem to be not sensitive to the definition of the stress tensor, at least 
qualitatively [TT1[S3]. We will discuss another possible definition of the stress tensor, later. 



2.3 Dynamics 

While we have specified the equilibrium probability distribution in the previous section, the 
dynamical properties such as the viscoelasticity cannot be determined unless the time evolution 
equations (rules) for the state variables are specified. In this section, we design time evolution 
equations which satisfy the detailed balance condition. (The detailed balance condition is required 
to be satisfied to reproduce the equilibrium thermodynamical properties correctly.) There are many 
possible time evolution equations (rules) which satisfy the detailed balance condition, and thus we 
cannot uniquely determine the dynamical model solely by the equilibrium probability distribution 
and the detailed balance condition. In this work, we propose a simple dynamics model which is 
similar to the single chain slip-spring models [10(111] and suitable for numerical simulations. 

First, for the time evolution of the bead position iZj^fe, we employ an overdamped Langevin 
type equation of motion. In absence of the external deformation field, the dynamic equation is 
given as follows. 

-^r- = dR-, ^^^^ 

Here C, is the friction coefficient of a bead and £,i.k{t) is the Gaussian noise obeying the fiuctuation 
dissipation relation of the second kind; 

(e:,fe(0>=0 (17) 

{£..Mmu{t')) = '^-^5.,,5ui5{t - t')l (18) 

where (■ ■ ■ ) represents the statistical average. For some analyses, it would be convenient to 
introduce the Fokkcr-Planck equation. The Fokker-Plank equation which corresponds to cq (|16p is 
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Here P{{Ri,k},{Sa}, Z;t) is the time dependent probability distribution and £fp is the Fokker- 
Planck operator. It is clear that eq ((T9)) satisfies the detailed balance condition and the steady 
state distribution coincides to the equilibrium distribution given by eq 

Second, we consider the reconstruction process of slip-springs. We assume that the reconstructions 
of slip-springs are independent of each other, and the positions of polymers and other slip-springs 
do not change during the reconstruction process. We write the construction rate of a new slip- 
spring as W+{Sz, Z\Z — 1) and the destruction rate of the /3-th slip-spring as W-{Z — 1|S'^, Z). 
We assume that a slip-spring is destroyed with a certain fixed probability when one of its ends is 
at the chain ends. The destruction rate can be written as 

W^{Z - l\Sp, Z) = ^ [55,,.,i + 5s,,,,N + <^s.,4,i + 5s,,,,n] £{li. Z) (20) 

Here, is the friction coefficient of a slip-spring and £{l3,Z) is the exchange operator which 
exchange the /3-th and Z-th slip-springs. (By operating £(/3, Z), and Sz are exchanged while 
the other slip-spring indices are unchanged. This operator is employed to ensure that the Z-th 
slip-spring is always destroyed.) 
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The slip-spring construction and destruction processes should be detailed-balanced. The detailed 
balance condition can be explicitly written as follows. 



W+iSz, Z\Z - l)P,qi{R,,k},{Sc.}, Z - 1) - l|5^,Z)Peq({H»,fc},{5a},^) (21) 

13 = 1 

From cqs ([20)) and pT|) . the construction rate is uniquely determined. By substituting eq (|20)) into 
eq (PT|) . we obtain the following explicit form for the construction rate. 
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As before, it would be convenient to introduce the dynamic equation for the time dependent 
probability distribution. The master equation for this reconstruction is written as 
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= C P 

where we have introduced the time evolution operator for the reconstruction process £ic- 

Finally, we consider the hopping of slip-springs along the chain. We assume that there is no 
interaction between slip-springs and each hopping event is statistically independent. The hopping 
process can be described by the change of connectivity index. For simplicity, we also assume that 
the change of connectivity index is restricted as ±1. (Namely, in our model, the hopping distance 
of slip-spring on a chain corresponds to the bead size). We consider the event where the (3-th 
slip-spring changes its connectivity from Sfj to 5^ . We describe the set of slip-spring indices after 
the hopping as {S^^}, for convenience. S'^ is defined as 



" \Sa (otherwise) 



(24) 



Let us indicate the transition rate of /3-th slip-spring from 5'^ to 5'^ as W{S'^\Si3) and its transition 
rate of the inverse process as VF(S'^|5'^). The detailed balance condition can be written as 



WjS'^lSp) ^ P,^{{R,^,},{S'J,Z) 
WiSp\S') P,q{{R,,k},{Sa},Z) 



exp 



(25) 



If we employ the Glauber type dynamics for the hopping process |56] . the transition rate which 



7 



satisfies eq (p5|) can be expressed as follows 
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The master equation for the hopping can be written as 
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Here we have introduced the time evolution operator for the hopping process £hop- 

Full dynamics of the system is then described by the Langevin equation given by eq (|16p . the 
reconstruction process (with the reconstruction rates (j20p and (|22p ). and the hopping process (with 
the hopping rate (|26p ). All of these processes satisfy the detailed balance condition, and thus it is 
clear that the equilibrium state is realized as characterized by the probability distribution ([5]) . The 
master equation of the system is expressed by combining the Fokker-Planck and master equations, 
dUl), and dMl)- 



dPi{R,,k}ASc.},Z;t) 
dt 

2.4 Relaxation modulus 



[£fP + -Crc + -Chop] 
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In this section, we will derive an explicit expression of the relaxation modulus tensor from the 
equilibrium distribution (eq ([H) and the master equation (eq ([Ml)) via the linear response theory 

m- 

We consider the system is subjected to the weak external deformation field, which is characterized 
by the time-dependent velocity gradient tensor K(t). Such a deformation gives an additional term 
to the master equation (P^)) . which will be treated as a perturbation in the foUowings. By adding 
the perturbation term, the master equation (j29p is modified as 



dP{{R,,k},{Sa},Z;t) 
dt 



= [Co + Ci{t)] P{{R,m}, {Sc.}, Z- t) 



where £0 and >Ci(t) are the equilibrium and perturbation time evolution operators. 
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(31) 
(32) 



As we mentioned, our dynamics model satisfies the detailed balance condition and thus the following 
relation holds for £0 and the equilibrium distribution (|5]) . 



UP,c{{R^..k}ASc.],Z)^Q 

Up to the first order in the perturbation, eq (j30p can be formally integrated as 



(33) 



P{{R^M}ASo.},Z-t) = P,^{{R,.k}ASo.},Z)+ f dt'e(*-*')^« [A(t')^cq({i?<.fe},{^4,^)] (34) 
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Then the ensemble average of the stress tensor at time t, cr{t), can be written as 
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Here we have defined the time shifted operator as a-{t) = e*^"o- {Cq is the adjoint operator for 
Co). This represents the stress tensor at time t after the reference time. Also, we have defined the 
virtual stress operator 6-^"' as 
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The virtual stress represents the stress generated by slip-springs and the repulsive interaction 
between beads. The relaxation modulus tensor G{t) (which is a fourth order tensor) can be 
defined for a small deformation as follows. 



o-(t) - (6-)cq = /" cW G{t - t') : K{t') 

J — OO 



By comparing eqs ([55]) and ([T7|) . we obtain 
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(38) 



Eq ([55]) is similar to the linear response formula obtained for single chain models [TTJ[3T] where 
the necessity of the virtual stress tensor is already known. However, the explicit form of the 
virtual stress tensor (p6|) differs from one for single chain models. In our model, the virtual stress 
tensor has the contribution from the repulsive interaction between beads. Physically this is natural 
because the repulsive interaction originates as the compensation of the attractive interaction by 
the slip-springs. 

In eq ((38l) . we assume that only the stress tensor a represents the stress tensor of the system 
(eq (|14p '). As wc mentioned, this assumption is based on the stress-optical law, which is empirically 
known to hold for various polymeric materials. However, from the view point of the virtual work 
method, it is also possible to employ a + o-^^^ (which is conjugate to the deformation) as the stress 
tensor of the system. If we employ the latter expression, the ensemble average in eq (j38p is replaced 
by {[a[t)+a^^\t)][& + 



Git) 



oq. Then we have the following formula. 
V 



knT 



{[a{t)+a^^\t)] [a + a^^)]) 



cq 



(39) 



(In the single chain slip-spring model, both of these two different expressions give qualitatively 
similar relaxation moduli. We expect that the situation is similar in our multi chain model.) 
Wc will compare simulation results for eqs ([M)) and (15^1) . later. 
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2.5 Numerical scheme 



In this section, we show a numerical scheme for simulations based on our multi-chain slip-spring 
model. We choose the bead size b, ksT, and C as the unit of length, energy, and friction (so that 
the unit of time is tq ~ CJj^ /ksT). In the foUowings we set 6=1, fc^T = 1, and C = 1- By using 
the operator-splitting method, the formal solution of eq (j29p can be approximated as 



P{{Ri,u)ASo.},Z-t + M) 



P({H,,fc},{5„},Z;t) 



(40) 



Eq (j40p corresponds to a numerical scheme with three substcps and the integration time step At. 
That is, the time evolution from time t to time t + At is simulated by performing the Langevin 
dynamics of the beads, the hopping dynamics of the slip-springs, and the reconstruction of the 
slip-springs. (These steps arc iterated sequentially.) 

For the integration of the Langevin equation, we employ the explicit Euler scheme. The 
Langevin equation for the chain dynamics can be discrctized as 



i2,,fe(t + At) « - At 



dT{{R,^k},{Sc.},Z) 



dR 



H,k 



(41) 



Here Wi^k is the Gaussian random number vector. 

The hopping dynamics of slip-spring is described by the change of slip-spring indices as mentioned 
above. When At is sufhciently small, the index of the a-th slip-spring, Sa,\ (A 2,4) is changed 
as Sa,x — > Sa,\ ± 1 by the following cumulative probability. 



At 



1 — tanh 



(A = 2, 4) 



(42) 



Here AJ^^i is the free energy difference given by 



2± 



2N. 



[{Rs 



2±1 



- RS^.3. 



4±1 



r- 



RSa..3:Sc,i) ] 

Rs - 



(43) 
(44) 



The reconstruction of the slip-springs is performed as follows. When an end of the a-th slip- 
spring is at a chain end, the slip-spring is destroyed with the following cumulative probability. 



At 

cT 



(45) 



When the slip-spring is destroyed, the index a for the other slip-springs is rearranged to realize 
a = 1,2,3,...,Z without vacant number. (This rearrangement is expressed by the exchange 
operator in eq ([20| .) After the attempts of destruction for all slip-springs, construction of a new 
slip-spring is attempted. This construction step is made by a Monte Carlo sampling scheme. A new 
slip-spring is virtually generated and its end is attached to a chain end. Another end is attached 
to one of the surrounding beads which is chosen randomly. Thus a new slip-spring index S is 
generated. This virtual slip-spring is accepted as a newly generated slip-spring with the following 
probability. 

3(i2si,S2 



exp 



-^53,54)^ 



2Ns 



(46) 



If accepted, we set Sz+i = S and increase Z ^ Z + 1. This Monte Carlo sampling is made for 



K = A^M^N^e" 



(47) 
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times on average, where the factor AM'^N is the total number of possible connections factorized by 
the number of ends for a slip-spring, 2, the total number of chain ends, 2M, and the total number 
of beads, MN . As a simple scheme, we assume that K only takes the floor or ceiling of {\.K\ 
or [/v ] ) . The probability that we have K trials at a certain construction step is given by 



P{K 
P{K 



[if] - K 

K - [K\ : 



(48) 
(49) 



This guarantees that the average sampling number becomes {K) = K, and the numbers of the 
slip-spring construction and destruction balance in equilibrium. 

The number of Monte Carlo sampling for construction can be significantly reduced if we exclude 
the constructions for very low acceptance probabilities. The acceptance probability decreases 
considerably if the stretch of a slip-spring becomes large. Thus, we limit the newly constructed 
slip-springs to be sampled only inside a certain cut-ofi^ size. We introduce a cut-off Tc = 
to restrict the sampling for construction of slip-spring with its length shorter than Tc- The cut-off 
parameter Co is determined to obey e~^°^'^ <^ 1 and practically this condition is satisfied if Co ^ 3. 
The average number of trials for every At is written as 



(50) 



Here Nq{Co) is the average number of beads located inside the distance CqNs/3 from the 
subjected chain end. No{Co) is written as 



J\r\<r,_ \fc=l 



N-1 



1 + E 



k=l 



4^ fC^oNsV^' 



Po 



N 
V 



exp 



47rr3 (M - 1)N 



V 



2k 



(51) 



Eq ((5T|) can be numerically evaluated if the parameters (such as Cq and Ns) are given, and as long 
as the parameters are unchanged during the simulation, K can be treated as constant. With this 
cut-off, the number of trials (for constant po) is 0{K) = 0{M) which is much smaller than that 
without cut-off, 0{K) = 0{AP). 



2.6 Comparison with earlier models 

In this section, to clarify the position of the proposed model, we compare our model with couples 
of similar models for entangled polymers. 

It has been rather established that Krcmcr-Grcst type coarse-grained molecular dynamics 
simulation [TH] is the standard way to reproduce polymer dynamics including entangled systems. 
In this approach, the multi-chain dynamics is solved with the excluded volume interaction that 
guarantees uncrossability between chains. On the other hand, in our model the interaction between 
chains is very soft and the entanglement effect is introduced a priori by the slip-springs. Since the 
relation between the contacts among chains in Krcmcr-Grcst simulations and the entanglement 
used in the entanglement-based models (such as our model) has not been clarified, our slip- 
spring reconstruction rules can not be related to the dynamics in Kremer-Grest simulations at 
this time being. For instance, it has been reported that in Krcmcr-Grcst simulations the long-lived 
contacts between two chains are constructed not only around the chain ends but also interior of the 
chain [58ll59] . But in the presented model (and the other entanglement-based models) assumes that 
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the entanglement reconstruction occurs only at the chain ends. It should be remarked, however, 
the resultant chain dynamics is similar to each other as shown later. 

Padding and Briels [33] have proposed a smart approach (referred as the TWENTANGLEMENT 
model) to deal with the uncrossability among chains without the excluded volume interaction. 
In the TWENTANGLEMENT model, a crossing event between bonds (which polymer chains 
consists of) is mathematically detected, and if the crossing occurs the force between segments is 
generated, to avoid chains cross freely. Due to the absence of the excluded volume interaction, 
their model is much coarse-grained than Kremer-Grest model and rather close to our model. On 
the other hand, the basic idea on the entanglement is the common for the Kremer-Grest and 
TWENTANGLEMENT models. Namely, both models do not require any artificial objects which 
represent the uncrossability (such as slip-springs in our model). We should also mention that there 
are no artificial attraction between chains in the TWENTANGLEMENT model, and thus the 
repulsive interaction is not required unlike our model. In these aspects, the TWENTANGLEMENT 
model is located in between the Kremer-Grest model and our model. 

Masubuchi et al |34j have developed another multi-chain model called the primitive chain 
network (PCN) model as mentioned in the introduction. The model presented in this work and 
the PCN model are similar to each other; in both models, phantom chains are connected to form 
a network, and the entanglement is mimicked by bundling of segments rather artificially. The 
reconstruction of entanglement is assumed to occur only at the chain ends as considered in the 
tube theory. One fundamental difference is the level of description. In the PCN model, only the 
number of Kuhn segments between entanglements is used and the position of each segment is not 
monitored. Thus, the PCN model cannot deal with the dynamics and structure in the time and 
length scales below those of entanglement. On the contrary, in the present model the dynamics 
of segments between entanglements is considered explicitly. This difference gives a difference 
in computational costs. The computational cost of the PCN is much smaller than one of the 
present model, as shown later. Another difference is the thermodynamic consistency, which is fully 
considered in the present model while not in the PCN model. The reconstruction of entanglements 
in the PCN model does not fulfill the detailed balance condition, for example [51]. Finally, we point 
that the strength of dynamical constraint is different; the slip-link employed in PCN corresponds 
to the limit of iVg = for the slip-spring of the present model. This difference may affect some 
properties such as the orientation tensor under deformations |43j . 

As mentioned in the introduction, the presented model is the many chain version of the 
Likhtman's single chain model [10] where several slip-springs are connected to a single Rouse 
chain. In the Likhtman's model, one end of a slip-spring can slide along the chain contour while 
another end is fixed in space. It the present model, on the other hand, both ends of a slip-spring 
can slide along chains. As we shall discuss later, this difference affects the effective strength of the 
dynamical constraint to the chain dynamics. Another difference between the two models is that 
the number of springs (standing for entanglements) is assumed to be constant in Likhtman's model 
while it fluctuates in the present model since it is controlled via the chemical potential. It is also 
pointed that there are a couple of differences in the sliding rule for slip-springs along the chain. As 
we mentioned, our slip-spring hops between segments while Likhtman's slip-spring actually slides 
on the bond between segments. (However, the difference between the hopping on beads and sliding 
along the bond between segments seems to have minor effects to dynamical properties, as judged 
from single chain simulation results |llj.) Furthermore, in Likhtman's model the slip-springs are 
allowed neither to change their order along the chain nor to overlap with each other, while in our 
model the motion of slip-spring is completely independent. 

It should be noted that there have been proposed several single chain slip-link models with 
the thermodynamic consistency. Schieber et al [211101111] have proposed such models where the 
state variables are chosen as the number of entanglement, the position of entanglement and the 
number of monomers between entanglements. This choice of state variables is similar to one of 
the PCN model, but fundamental difference is that these models have well-defined free energy free 
energy of the system. The dynamic equations or transition rates are derived from the free energy 
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and the detailed balance condition. The chemical potential controls the fluctuation of number 
of entanglements and this strategy is employed in our model. Due to the nature of the single 
chain model, it is intrinsically difficult to deal with the effect of surrounding chains. (Although 
several attempts have been made to overcome this difficulty |62[63j , some additional and non-trivial 
assumptions are required to mimic multi chain effects.) 

2.7 Simulations 

Monodisperse linear polymers were examined where bead number per chain was varied from 4 
to 64. The total number of beads in the system, M N, is fixed to be constant so that the bead 
number density was constant at po = 4. Periodic boundary condition was utilized with the box 
dimension at 8"^. The spring strength parameter Ng for slip-springs was set as Ns ~ 0.5. The 
number density of slip-springs was chosen as (/> = 0.5. Conceptually, this slip-spring density give a 
certain plateau modulus but we have not yet obtained the relation between the plateau modulus 
and these parameters. The cut-off parameter and the corresponding cut-off distance were Cq = 10 
and = 1.29, respectively. The friction coefficient for the slip-spring was set as Cs = C- For the 
numerical calculation At was chosen as 0.01, after we checked reasonable numerical convergence 
for At < 0.02. 

3 Results and Discussion 

3.1 Statics 

In this subsection static properties of the system is examined to show the consistency between 
simulation data and the equilibrium distribution function given by eq Figure [2] shows the bond 
number {N — 1) dependence of squared end-to-end distance to report that the scaling obeys 
the Gaussian chain statistics. The internal chain structure is examined via the internal distance 
factor d{s) defined as d{s) = {{Ri^k+s — Ri,k)'^)/s, and shown in Figure [H d{s) is reasonably 
close to unity independently of s and N as expected for the Gaussian chain statistics. These 
results demonstrate that the attractive force induced by the slip-springs is correctly compensated 
by the repulsive interaction. A similar attempt has been made for PCN model where the soft 
core repulsive interaction was introduced between slip- links |64| , but a precise control of the chain 
dimension was difficult due to the lack of free energy expression for PCN model. (Indeed, it 
seems practically impossible to introduce the repulsive potential which exactly cancels the artificial 
attractive interaction in slip-link models such as the PCN model [53].) 

Figure m shows the distribution of slip-spring number per chain Zc, P{Zc). The grand canonical 
type treatment of the slip-springs for single chain models predict the Poisson distribution for the 
number of slip-springs on a chain [TT1[54]. Indeed, the results shown by symbols reasonably close 
to the Poisson distribution drawn by solid curves where the average value of Zc is given as 2,4,8 
and 16 for the examined chains with A^ = 8, 16, 32 and 64, for the simulation parameter A'e = 4 
(calculated from the slip-spring density and the bead density as A^e ~ 2(/)po)- 

3.2 Dynamics 

Figure [5] shows the mean square displacement of the central bead in the chain, g^^'^{t). To see the 
effect of the entanglement clearly, the data divided by the Rouse behavior in internal time range 
g™^'^(t)t~^/^ is also shown. From these plots, it is apparent that the entanglement effect appears 
after t ^ 10 where the negative slope starts in (7™"^(<)i^^/^. In this respect, the short chain with 
A^ = 8 does not show the entangled behavior, in spite of the fact that there exist two slip-springs 
per chain on average as shown in FigU) Figure IH] shows the diffusion coefficient D against the bead 
number A^. To see the scaling behavior clearly, DN^ is also shown. The Rouse behavior D cx A^~^ 
is observed below N = 8, which is consistent with what observed for g™^'^{t). For longer chains. 
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the N dependence is not that strong and actually the power law exponent does not largely deviate 
from —2 even for N = 64, which is somewhat weaker than experimental results for well entangled 
systems j65j . In comparison with the literature data (shown by cross in the right panel) for single 
chain slip-spring model |10j with similar parameter setting (except the difference in C.s that is 1 
for our simulation while it is 0.1 for the single chain data), it is suggested that the dynamical 
constraint in the present model is weaker than that in the single chain model. We expect that this 
is due to the difference in the hopping (sliding) kinetics of slip-springs. In our model the both ends 
of slip-springs are mobile while in the single chain models one end is anchored, as we mentioned in 
Section 12.61 This difference will result in a weaker constraint to the chain in our model. Indeed, 
similar behavior is observed in the N dependence of the longest relaxation time for the end-to-end 
vector, Tmax, showu in Fig[71 Especially, Tmax.^"^ indicates that the large N chains examined here 
are still in the transitional region between unentangled and entangled regimes. 

Figure |8] shows the relaxation modulus G{t) calculated by the linear response theory shown 
in Sec. 12.41 Here, G{t) is normalized by the unit modulus defined as Go = pRT/Mo [66]. As in 
the case of the mean square displacement, we also show the data divided by the Rouse behavior, 
G{t)t^^'^. As we have observed for g™"^(i), the entanglement effect (the increase 
is not observed in short chains. On the other hand, for longer chains G{t) apparently deviates from 
the Rouse relaxation to show the plateau, as expected. We note that the unit modulus Go is not 
equal to the plateau modulus Gat and actually Gat is much smaller than Gq. (This is because Gq 
is a characteristic modulus for the bead scale, and it reflects the short time relaxation modes.) For 
the single chain model it is reported that Gn ~ O.lGo [IH]- Our results for long chains are similar 
to this relation. Of course, the relation between Gn and Go depends on the model and various 
parameters. The direct comparison of our result with the single chain simulation results seems to 
be difficult. 

As we mentioned, there are two possible expressions for the stress tensor of the system. To 
investigate how the expression of the stress tensor affects the rheological properties, here we 
compare the shear relaxation moduli calculated with two different expressions. Figure IHl shows 
the shear relaxation moduli for = 8, 16, 32 and 64 calculated by eqs ((38|) (symbol) and (p9)) 
(solid curve). Although the shear relaxation moduli data by eqs (p8| and ([39]) are quantitatively 
different, they are qualitatively similar, as reported in the single chain model. |11| 

Figure ITOl shows a comparison with the Kremer-Grest simulations |66] on G{t). Here, G{t) 
of our simulation is calculated by eq but similar fitting can be realized for G{t) obtained 

by cq (p9|) as well, if the parameters are adequately tuned (not shown here). There are small 
discrepancies in the short time region due to the difference in the level of coarse-graining, i.e., the 
number of Rouse modes. This discrepancy may be eliminated if A'^e value is increased as reported 
in the single chain model. For this fitting, we choose the parameters as follows. 

Afo = 3.1, TO-50TKG, Go = 0.2Go,KG (52) 

Here, Mq is the number of beads (molecular weight) of the Kremer-Grest model which corresponds 
to one bead in our model, tkg is the unit time in the Kremer-Grest simulation (the standard 
time scale for particles with the Lcnnard- Jones potential), and Go,kg is the unit modulus for the 
Kremer-Grest model |66| . This fitting demonstrates that the dynamical constraint given by the 
slip-spring is relatively weak, and this result is consistent with results obtained by D and Tmax- For 
instance, it has been reported that the entanglement number for the Kremer-Grest chain with 200 
beads is around 5 |66| . On the other hand, the corresponding chain in our model is A^ = 64 that 
has 16 slip-springs per chain on average, as shown in Fig |4l Since the G{t) data compared here 
are in the transitional regime and the contribution from the Rouse relaxation is relatively large, 
the fitting results (cq ([52]) ) may be different for well-entangled systems. Nevertheless, the fact 
that Ne is smaller than the effective entanglement bead number is informative. For this specific 
case, the effective entanglement bead number in our model is roughly estimated as 3A^e- A similar 
result have been reported for the PCN model [Sj, and such a discrepancy is expected to be the 
nature of randomly connected multi chain network. From the similarity between our model and 
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the PCN model, where the entanglements are not fixed in space, we consider the obtained results 
are reasonable. 

In order to compare the computation cost (efRciency) of our model, we report CPU times for 
several different simulation models on a PC (Intel Xeon X5570 2.93GHz, Cent OS 4). We chose a 
melt of Kremer-Grest chains with bead number N = 200 as a reference, and made calculations by 
the present model and by the PCN model for similar systems. Each simulation was performed for 
time steps comparable to the longest relaxation time. All the simulations were performed by serial 
simulation codes. The Kremer-Grest simulation was made by COGNAC ver 7.1.1 [67]. The other 
simulations were performed by our home-made simulation codes. For a Kremer-Grest simulation 
with the chain number M = 6, the required time steps is about 100, 000 and the CPU time was 60 
hours. For our model, a simulation with N = 64, M = 32, Ng = 4 and (s = 1-0, required the 3000 
time steps and the CPU time was 30 minutes. For a PCN simulation with the average segment 
number of segments per chain Z = 15 and the number of chains M = 341 required 200 time steps 
for the relaxation and the CPU time was 2 minutes. From these results, we can conclude that the 
computational cost (efficiency) of our model is in between the Kremer-Grest model and the PCN 
model. Judging from the coarse-graining level of these models, this result is reasonable. 

4 Conclusion 

In this study, we proposed the multi-chain slip-spring model where the bead-spring chains are 
dispersed in space and connected by slip-springs. The entanglement effect is mimicked by the 
slip-springs and not by the hard-core (excluded volume) interaction between beads. Our model is 
located in the niche of the modeling of entangled polymer dynamics between conventional multi- 
chain simulations and single chain models. In our model, the set of state variables are the position 
of beads and the connectivity (bead indices) of the slip-springs. Differently from the primitive chain 
network model (that is the other class of multi-chain model), the total free energy of the system is 
well-defined, and kinetic equations are designed based on the free energy and the detailed balance 
condition. The free energy includes the repulsive interaction between beads, which compensate 
the attractive interaction artificially generated by the slip-springs. The explicit expression of linear 
relaxation modulus was also derived by the linear response theory. A possible numerical scheme 
was proposed and simulations reproduced expected bead number dependence in transitional regime 
between Rouse and entangled dynamics for the chain structure, the central bead diffusion and the 
linear relaxation modulus. 

In this model there are a couple of parameters {Ng, Cs, Po and </>) of which effect on chain 
dynamics is not well understood. It would be an interesting work to explore the effects of these 
parameters by theoretical methods as well as simulations. We are performing simulations to scan 
these parameters and compare the resulting chain dynamics with the other models and experiments. 
Another interesting topic is on the extension of this model by using its multi-chain nature, as 
reported for PCN model, to branch polymers, blends and copolymers, etc. The results for these 
attempts will be reported elsewhere. 
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1 Schematic representation of the model Bead-spring chains are dispersed in space 
and connected via shp-springs with each other to form network 

2 Squared end-to-end distance plotted against bond number 

3 Internal distance factor d{s) for = 8, 16, 32 and 64 indicated by square, triangle, 
cross and circle, respectively ^ 

4 Distribution of slip-spring number per chain Zc for N — 8, 16, 32 and 64 indicated 
by square, triangle, cross and circle, respectively. Poisson distributions with the 



average values of 2, 4, 8 and 16 are shown by solid curves |21 

5 Mean square displacement of the central bead for N = 8, 16, 32 and 64 from left to 



right |22 

6 Diffusion coefficient D plotted against bead number N. Right panel shows comparison 
with the single chain model |10j indicated by cross 

7 Relaxation time for end-to-end vector Tmax plotted against bead number N 

8 Linear relaxation modulus for TV = 8, 16, 32 and 64 from left to right 

9 Linear relaxation moduli for = 8, 16, 32 and 64, calculated with two different 
expressions of the stress tensor (eqs ([38|) shown by symbol and (|39| shown by solid 



23 
23 
24 



line) 125 

10 Linear relaxation modulus in comparison with the Kremer-Grest simulations. For 
the Kremer-Grest simulations the bead numbers are 50, 100 and 200 and the data 
taken from Ref. [B^ are shown by filled symbols according to left and bottom axes. 
For the present model the bead numbers are 16, 32 and 64 from left to right and 



the data are shown by unfilled symbols according to right and top axes [25 
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Figure 1: T. Uneyama and Y. Masubuchi to J. Chem. Phys 
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Figure 2: T. Uneyama and Y. Masubuchi to J. Chem. Phys 
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Figure 3: T. Uneyama and Y. Masubuchi to J. Chem. Phys 
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Figure 5: T. Uneyama and Y. Masubuchi to J. Chem. Phys 
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Figure 6: T. Uneyama and Y. Masubuchi to J. Chem. Phys 
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Figure 7: T. Uneyama and Y. Masubuchi to J. Chem. Phys 
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Figure 8: T. Uneyama and Y. Masubuchi to J. Chem. Phys 
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Figure 9: T. Uneyama and Y. Masubuchi to J. Chem. Phys 



t/To 

12 3 4 

10 10 10 10 




10^ 10^ 10^ 10^ 10^ 

Figure 10: T. Uneyama and Y. Masubuchi to J. Chem. Phys 
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